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Abstract 

Background: Waddington's epigenetic landscape is an intuitive metaphor for the developmental and evolutionary 
potential of biological regulatory processes. It emphasises time-dependence and transient behaviour. Nowadays, we 
can derive this landscape by modelling a specific regulatory network as a dynamical system and calculating its 
so-called potential surface. In this sense, potential surfaces are the mathematical equivalent of the Waddingtonian 
landscape metaphor. In order to fully capture the time-dependent (non-autonomous) transient behaviour of 
biological processes, we must be able to characterise potential landscapes and how they change over time. However, 
currently available mathematical tools focus on the asymptotic (steady-state) behaviour of autonomous dynamical 
systems, which restricts how biological systems are studied. 

Results: We present a pragmatic first step towards a methodology for dealing with transient behaviours in 
non-autonomous systems. We propose a classification scheme for different kinds of such dynamics based on the 
simulation of a simple genetic toggle-switch model with time-variable parameters. For this low-dimensional system, 
we can calculate and explicitly visualise numerical approximations to the potential landscape. Focussing on transient 
dynamics in non-autonomous systems reveals a range of interesting and biologically relevant behaviours that would 
be missed in steady-state analyses of autonomous systems. Our simulation-based approach allows us to identify four 
qualitatively different kinds of dynamics: transitions, pursuits, and two kinds of captures. We describe these in detail, 
and illustrate the usefulness of our classification scheme by providing a number of examples that demonstrate how it 
can be employed to gain specific mechanistic insights into the dynamics of gene regulation. 

Conclusions: The practical aim of our proposed classification scheme is to make the analysis of explicitly 
time-dependent transient behaviour tractable, and to encourage the wider use of non-autonomous models in 
systems biology. Our method is applicable to a large class of biological processes. 



Dynamical systems theory, Differential equations, Non- 
autonomous systems, Transient dynamics, Genetic toggle 
switch, Phase space analysis, (Quasi-)potential landscape 

In the study of development we are interested not only in 
the final state to which the system arrives, but also in the 
course by which it gets there. 

Conrad Hal Waddington, "The Strategy of the Genes" 1957 
Background 

Development in wild-type organisms is robust to genetic 
and environmental variations. Conrad Hal Waddington 
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introduced the notion of canalisation to describe how 
developmental processes resist perturbations during 
embryogenesis [1-3]. The canalised nature of develop- 
ment explains, he argued, why most phenotypes are dis- 
crete and distinct. To illustrate these ideas, he developed 
the epigenetic landscape, one of his most well-known 
concepts [3,4]. 

In Waddington's epigenetic landscape, the current state 
of a developing system is indicated by a ball on an 
undulated surface (Figure 1A, top panel) [3,5,6]. The 
topography of this landscape determines the develop- 
mental potential or repertoire of the system. The top- 
most edge of the surface shown in Figure 1A (top panel) 
represents the initial state of the system given by, for 
example, a particular set of initial protein concentra- 
tions in a cell. Valleys in the landscape symbolise the 
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Figure 1 Waddington's epigenetic landscape and potential surfaces. (A) Two different views of Waddington's epigenetic landscape taken 
from "The Strategy of the Genes" published in 1 957 [3]. The top panel shows a top view of the landscape. The path that the ball will follow 
represents the developmental trajectory (chreode) of a given system. Valleys indicate alternative differentiation pathways, branch points imply 
developmental decisions. The bottom panel shows the view from below the landscape. It illustrates how genes remodel the surface by pulling on it 
through ropes. Waddington used this sketch to show how the landscape's topography changes during development and evolution. (B) Panel 1 : 
Diagrammatic representation of the toggle switch network used in the simulations. Activating interactions are indicated by arrows, repressing ones 
by T-bar connectors. See Model and methods section for detailed parameter descriptions. Panel 2: Mathematical formulation of the toggle switch 
model, x and y indicate concentrations of the protein products of genes X and Y. Ordinary differential equations define the rate of change in protein 
concentrations and ^). Sigmoid functions with fixed Hill coefficients of 4 are used to represent auto-activation and mutual repression. Decay and 
external activation are taken to be linear. Parameters as in Panel 1 . Panel 3: Phase portrait for a constant set of parameter values of the toggle switch 
model in the bistable regime. X- and /-axes represent protein concentrations x and y. We use this example to illustrate relevant features of phase 
space: arrows indicate flow, blue points mark the position of stable steady states (attractors), the red point shows an unstable steady state (saddle) 
lying on the separatrix that divides the two basins of attraction (grey line). See the main text for detailed descriptions of the highlighted features. 
Panel 4: Quasi-potential landscape associated to the the phase portrait shown in Panel 3. The steepness of the quasi-potential surface correlates with 
the flow at each corresponding point on the phase portrait. Attractors, saddle, and separatrix are indicated as in Panel 3. See main text for details. 



various differentiation pathways that are available. The 
landscapes topography — together with the initial state — 
determine a developmental trajectory that follows a par- 
ticular valley. The structure of the landscape is such that, 
if the system is slightly perturbed, the sloping valley walls 
will cause it to correct and readjust its trajectory. This 
behavour is called 'homeorhesis' — the maintenance of a 
dynamic trajectory — in analogy to the more static con- 
cept of homeostasis— the maintenance of a (steady) state 
of the system [3]. The wider and deeper a valley is, the 
more canalised the developmental trajectory. Waddington 
named such canalised trajectories chreodes'. 

An additional feature of Waddington's landscape is cru- 
cial in our context: it is the fact that its surface is not 
necessarily fixed, but can change over time due to the in- 
fluence of genetic or environmental signals. Waddington 
represented this by pegs connected to the underside of the 
landscape by ropes (Figure 1A, lower panel) [3]. As the 



genetic or environmental context of the system changes, 
these ropes pull and stretch the surface, thereby changing 
its topography and hence its developmental potential. 

Waddington's epigenetic landscape was conceived as a 
conceptual tool to illustrate the nature of developmen- 
tal robustness and its effect on evolutionary dynamics 
[1,3-5,7]. As such, it remained at a rather metaphorical 
level of explanation, since complex non-linear biological 
processes are hard to formulate and analyse [5,6,8-12]. 
This is still the predominant way in which it is used in 
various reviews on contemporary stem cell research (see 
for example, [13-17]). However, in order to understand 
the precise nature of specific developmental trajectories 
or chreodes in real systems, we have to take Waddington's 
landscape a step further: we have to calculate it based on 
experimental evidence, and use it to characterise the tran- 
sient behaviours that govern the observed developmental 
dynamics [18-25]. 
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The increasing availability of quantitative gene expres- 
sion data renders this approach feasible. However, before 
we can successfully apply it, we also require new concep- 
tual and mathematical tools to deal with the analysis of 
data-driven models that are formulated in terms of non- 
autonomous (L e. explicitly time-dependent) dynamical 
systems. Explicit time-dependence is necessary to repro- 
duce the changing topography of Waddingtons landscape. 
However, such systems are difficult to study in a rigorous 
mathematical manner, and few analysis tools exist at this 
point. In this work, we address this challenge by proposing 
a classification scheme for transient dynamic behaviours 
observed in a non-autonomous version of a simple gene 
network model. This scheme is meant to provide the foun- 
dation for the analysis of more complex time-dependent 
models that reproduce the dynamics of specific, experi- 
mentally tractable, biological regulatory systems. 

Dynamical systems 

In this study, we focus on dynamical systems formulated 
in terms of ordinary differential equations, and illustrate 
how such models can help explain the function and poten- 
tial of developmental gene regulatory networks in terms 
of their dynamical repertoires, that is, in terms of the set 
of behaviours that can be implemented by the system. 
A systems behaviour is defined by its trajectories, which 
represent the change of the state of the system— e.g. 
consisting of a set of transcription factor concentrations- 
over time. The shape of these trajectories depends on the 
structure or organisation of the underlying regulatory net- 
work (see, for example, [3,22-24,26-31]). It is possible to 
gain a general qualitative understanding of what the sys- 
tems trajectories can and cannot achieve without the need 
to solve the dynamical system analytically [32,33]. Most 
of this information comes from the geometrical analysis 
of phase space, i. e. an analysis of the number, nature and 
relative arrangement of the steady states of the system. 

The phase space of a dynamical system is an abstract 
space, in which each dimension represents the value of 
a specific state variable. Here we use a well established 
double-repressive feedback loop model, known as the tog- 
gle switch ([23], and references therein), to study transient 
dynamics in a two-dimensional, time-dependent gene reg- 
ulatory network (see Figure IB, panels 1 and 2, as well as 
'Model and methods' below, for the full formulation). In 
this case, the state variables represent the concentrations 
of the transcription factors that constitute the network 
(denoted by X and Y; Figure IB, panel 1). The graphical 
representation of phase space is called the phase portrait. 
It shows the rate of change of the system at any given 
state. This is known as the flow of the system. The flow 
of the toggle switch model is indicated by arrows of a 
given length and direction in the phase portrait shown in 
Figure IB, panel 3. If we follow the flow from all possible 



initial states, we obtain the totality of possible dynamic 
trajectories. 

It is evident from the inspection of the flow in Figure IB 
(panel 3) that trajectories tend to converge to specific 
points in phase space: the steady states of the system. 
There are different kinds of steady states, those that are 
stable, and those that are unstable. The most simple sta- 
ble steady state is an attractor point [32,33]. Attractors, 
as their name implies, draw trajectories towards them. 
Furthermore, they have the special property that once a 
trajectory has reached an attractor, it will return to it if the 
system is slightly perturbed. An example of an unstable 
steady state is a saddle point (Figure IB, panel 3). Saddle 
points attract trajectories from some directions, but repel 
them in others. Usually, the system will move away from 
a saddle upon perturbation, towards the nearest attractor. 
The repelling trajectory follows a structure in phase space 
called an unstable manifold. This manifold is defined by 
the the path that links a saddle with an attractor point. 
Note that unstable manifolds correspond to chreodes at 
the bottom of valleys in Waddingtons landscape. 

In our phase portraits, we plot steady states as points 
coloured according to their stability: attractors in blue, 
and saddle points in red (Figure IB, panel 3). The region 
of phase space around an attractor, from which all trajec- 
tories converge towards it, is called its basin of attraction. 
Curves known as separatrices set apart the different basins 
and their attractors (in the case of Figure IB, panel 3, 
the separatrix is a straight line indicated in grey). Saddle 
points are always located on separatrices. 

Attractors and saddles, with their associated basins and 
separatrices can be created or annihilated through the 
process of bifurcation [32,34]. Bifurcations represent sud- 
den qualitative changes in the structure of the phase 
portrait caused by small changes in the values of a given 
set of control parameters. 

To date, the best example illustrating the importance 
and usefulness of geometric approaches for understand- 
ing the dynamics and function of specific, experimentally 
tractable, developmental regulatory systems comes from 
an analysis of the gap gene regulatory network involved 
in pattern formation during early embryogenesis of the 
vinegar fly Drosophila melanogaster *. Manu and colleagues 
[35] identified features of phase space responsible for pat- 
terning and canalisation of spatio-temporal gene expres- 
sion in the Drosophila blastoderm embryo. The analysis 
is based on low-dimensional projections of phase space 
to study the geometric arrangement of attractors in the 
four-dimensional system representing the change in pro- 
tein concentration for four transcription factors encoded 
by the gap genes hunchback (hb), Kruppel (Kr), knirps 
(kni), and giant (gt). Gene expression domain boundaries 
that remain at a constant position over time could be 
attributed to movements of attractors in phase space, or 
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nuclei switching between attractor basins depending on 
their position in the embryo. In contrast, more posterior 
expression boundaries, which keep shifting position over 
time, were associated with transient behaviours along a 
canalising unstable manifold, forming the equivalent of a 
valley in Waddingtons landscape [35]. 

Potential landscapes 

Phase portraits of systems with two state variables can be 
visualised in terms of their associated potential landscape 
[32] (Figure IB, panel 4). In this representation, the steep- 
ness of the potential landscape corresponds to the flow 
of the system (compare Figure IB, panels 3 and 4). Tra- 
jectories travel downhill towards the attractor points. The 
topography of the potential landscape is therefore a direct 
result of the topology (number and nature of steady states) 
and the geometry (relative positions of the steady states 
and size of the flow) of the underlying phase portrait. 

Potential surfaces can be thought of as mathemati- 
cal representations of Waddingtons epigenetic landscape 
[18,20,24,25,36-40]. Attractors represent differentiated 
states (in agreement with earlier postulates, [1,13,14,41-46]), 
separatrices form the ridges between the valleys, which 
are formed by unstable manifolds representing their asso- 
ciated chreodes (see Figure IB, panels 3 and 4), and the 
flow is represented by the steepness of the slopes on the 
landscape. Canalisation is explained by the depth and 
width of each valley in the potential landscape. This pro- 
vides a way of probing and understanding the features 
that confer robustness to the system (see, for exam- 
ple, [18,19,21,24,45,47,48]). Branching valleys can arise 
through particular types of bifurcation events. In particu- 
lar, new attractor states can branch off from existing ones 
during development or stem cell differentiation through 
(supracritical) pitchfork bifurcations [19,20,25,32,47,49-51]. 

Due to this analogy to Waddingtons epigenetic land- 
scape, potential landscapes are becoming increasingly 
popular as explanatory tools in fields such as evo- 
devo, developmental biology, and especially in stem cell 
research. In the case of stem cells, the positioning of the 
valleys in the landscape relative to each other explains 
which differentiation pathways can be reached by a given 
stem cell, and which cell fates can (or cannot) be trans- 
differentiated into each other [13-16,36,37,40,52-54]. 
This illustrates how specific biological meaning can be 
gained from studying the topography of potential land- 
scapes, which ultimately, is nothing more than a visu- 
ally accessible way of studying the underlying phase 
portrait. 

Waddington meant his epigenetic landscape as a "dia- 
grammatic representation" of development and warned 
explicitly against interpretations that were too rigorous 
or literal [2,3]. Detailed topographical interpretations of 
Waddingtons landscape may be quite inaccurate and even 



misleading. Ferrell [25] points out that Waddingtons land- 
scape, where valleys progressively split into an increasing 
number of branches, does not help explain many realistic 
cell differentiation processes. In particular, many induc- 
tive processes in development (e. g. vulval induction in the 
roundworm Caenorhabditis elegans, or mesoderm induc- 
tion in vertebrates) involve saddle-node rather than pitch- 
fork bifurcations, which correspond to the disappearance 
of valleys rather than to their creation by branching. The 
system shifts to a new attractor only once the old one 
has vanished. This is why it is important to move from 
metaphorical uses of Waddingtons epigenetic landscape 
to accurately calculated potential surfaces whenever pos- 
sible, or more importantly, to the detailed analysis of the 
underlying phase portraits, which is where the dynamics 
of specific regulatory networks are determined. 

Potential landscapes can only be calculated and visu- 
alised explicitly for a restricted range of dynamical sys- 
tems, belonging to the class of gradient systems [32]. Note 
that the notion of a gradient system is defined math- 
ematically by the absence of limit cycles or any other 
complex attractor structures in their phase portraits, and 
has nothing to do with biochemical or other biological 
gradients (see Model and methods for a detailed explana- 
tion). In cases where we do not know whether the system 
under study is a gradient system or not, we can approxi- 
mate the actual potential using various numerical meth- 
ods [19,20,38,39,48,55]. The resulting approximations are 
called quasi-potential landscapes. 

However, even such quasi-potential landscapes can only 
be visualised directly when the number of state vari- 
ables of the system does not exceed two. This is not 
true for most biologically realistic systems. Nevertheless, 
(quasi-)potential landscapes are still useful as conceptual 
tools for the analysis of higher-dimensional regulatory 
networks. In some cases, it is possible to reduce a high- 
dimensional system to a lower-dimensional one (see, for 
example, [56-58]). But even if this is not the case, the 
concept of the potential landscape provides two advan- 
tages. First, as discussed above, quasi-potential surfaces 
link Waddingtons intuitively accessible concept of the epi- 
genetic landscape to the biological interpretation of (high- 
dimensional) phase space analysis. And second, potential 
surfaces are useful as visual guides to diagnose features of 
the underlying phase portrait of the system that are char- 
acteristic of specific dynamic behaviours of the regulatory 
network under study. 

The importance of transient dynamics 

Many models of biological systems are formulated with 
the assumption that the relevant dynamics occur near or 
at a steady state. For instance, Thorns pioneering sys- 
tematic and rigorous analysis of morphogenesis in terms 
of catastrophe theory explicitly and strongly relies on 
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this assumption [8]. Similarly, work on robustness and 
evolvability, using ensembles of simulated networks, stan- 
dardly assumes that the steady state pattern produced by 
a model can be taken as a satisfactory and realistic rep- 
resentation of the phenotype of the system (e. g. [59-62]). 
Furthermore, stem cell models that make use of potential 
landscapes are analysed with a strong focus on how their 
attractors govern system behaviour [19,20,47,48,56,58]. 

In some cases, the steady state assumption makes obvi- 
ous sense based on biological reasoning. One exam- 
ple is the analysis of the segment polarity network in 
Drosophila melanogaster, which amplifies and maintains 
a periodic input, and thus performs an intrinsically sta- 
bilising patterning function [63-65]. In most cases, how- 
ever, biological pattern formation is highly dynamic and 
far from equilibrium, and the steady state assumption 
is justified based on methodological, rather than biolog- 
ical considerations. Focussing on asymptotic behaviour 
at or near steady state greatly simplifies the analysis of 
the system. First, it discretises and reduces phase (and 
hence phenotype) space into a small number of pos- 
sible states— represented by the systems attractors (see 
Figure 1C, panel 3). Second, it enables the powerful toolkit 
of linear stability analysis to be employed to examine the 
characteristic properties of system states [32,33]. 

However, there are both theoretical and practical rea- 
sons indicating that steady state analysis misses essential 
and biologically relevant systems behaviours. One line of 
theoretical reasoning is provided by Waddington himself, 
who reminds us that "[i]n the study of development we 
are interested not only in the final state to which the sys- 
tem arrives, but also in the course by which it gets there". 
[3]. His concept of homeorhesis and his representation of 
developmental trajectories or chreodes as descending val- 
leys in the epigenetic landscape places the focus explicitly 
on the transient dynamics of cells on their way to their 
final, differentiated state. In the same spirit, other authors 
have suggested that phenotypes should be defined over 
developmental trajectories, rather than representing some 
sort of 'final' outcome or end state' of the system [22,66]. 
To decide what a final pattern is, and when exactly it 
occurs, is always arbitrary to some degree, while transient 
features (such as intermediate stages and the timing of 
their transitions) are clearly important when considering 
the function and dynamics of developmental processes. 

There are also practical reasons to consider transient 
dynamics explicitly. For developmental processes that 
consist of continuous transitions between patterns rather 
than the production of a final output, it is impossible 
to decide a priori whether the system is representing a 
non-autonomous succession of steady states (see below), 
or whether its behaviour is truly transient (i. e. far from 
steady state). In the case of gap domain shifts, we have 
evidence for the latter [35], although there is no reason 



to assume that the two situations need to be mutually 
exclusive. The gap gene model analysed by Manu and 
colleagues [35] exhibits boundary shifts that are caused 
by trajectories following a canalising unstable manifold. 
Assuming steady state dynamics would collapse the trajec- 
tories representing these shifts into a single attractor point 
at the final configuration of gene expression. Any such 
analysis would miss the relevant underlying features of 
phase space (the transient manifold), and therefore fail to 
provide a proper characterisation and explanation for the 
observed gene expression dynamics (the shift in domain 
position over time). 

Other examples of developmental processes, where 
transient dynamics are clearly important are dorso-ventral 
patterning of the vertebrate neural tube, which involves 
boundary shifts strikingly similar to those of the gap 
domains [67,68], and vertebrate somitogenesis or short- 
germband segmentation in arthropods where transient 
travelling waves of gene expression are an essential com- 
ponent of the underlying clock-and-wavefront patterning 
mechanism [69-75]. Incidentally, similar considerations 
can be made for models dealing with ecological networks, 
and several examples exist in the literature that consider 
transient dynamics explicitly (see models of coupled oscil- 
lating population dynamics between species [76,77], or 
[78-80]). 

Non-autonomy: explicit time dependence 

Another important aspect of biological regulatory pro- 
cesses, which receives surprisingly little attention, is the 
explicit time dependence of these systems. As soon as 
we consider cellular dynamics, development, or evolution 
over a large-enough time span, the organisation of the 
underlying regulatory system starts to change. This affects 
the parameters — not just the state variables — of the sys- 
tem. Such explicitly time-dependent dynamical systems 
are called non-autonomous [32,81,82]. Time-dependent 
signalling cues and environmental conditions have long 
been known to shape many processes in the fields of evo- 
lutionary and developmental biology. Obvious examples 
of such phenomena are inductive processes or external 
(e.g. seasonal) cues that are essential to trigger many 
developmental pathways (as described in standard text- 
books such as [83,84]), or evolutionary dynamics driven 
by changing environmental conditions (examples, based 
on the simulation of gene regulatory network models, can 
be found in [85-88]). 

Still, it is rare to find studies based on explicitly non- 
autonomous models in the literature, and most authors 
avoid the challenge of dealing with dynamical systems 
where the parameters representing external cues are time- 
dependent. This is the case in the study of gap domain 
shifts by Manu and colleagues [35,89], where maternal 
morphogen gradients providing regulatory input to the 
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system were assumed to reach steady state before gap gene 
boundary positioning was analysed. Such simplifications 
can be risky, especially when describing biological phe- 
nomena where the time scales of change in parameters 
and state variables are of similar order. In such cases, time 
scales should not be separated, nor quasi-steady states 
considered since it is easy for dynamical properties and 
behaviours of the system to be missed or misinterpreted 
under these conditions. 

Recently there have been some attempts at including 
non-autonomy in biological models. Corson and Siggia 
[90], for example, offer an explanation for vulva devel- 
opment in C. elegans which explicitly considers temporal 
parameter changes due to inductive signalling cues in 
their model. Their model considers differentiation of cells 
into three differentiated states depending on inputs from 
two signalling pathways. Signalling inputs are encoded 
in the model by altering values of system parameters, 
which change and distort the geometry of the phase 
portrait by displacing separatrices from their original 
position (see Figures three-nine in [90]). When a cell 
receives a signal, its developmental trajectory comes to lie 
within another basin of attraction, inducing an alterna- 
tive cell fate to that which would have been reached in the 
absence of the signal. This study illustrates the importance 
of non-autonomous dynamics in development. However, 
it remains somewhat limited in its implementation of 
explicit time-dependence. Although parameter change 
is included in the model, signalling occurs before cells 
embark on their developmental paths, and trajectories 
develop purely autonomously thereafter. Therefore, this 
approach does not fully capture the transient dynam- 
ics of cell differentiation. In other words, although their 
model is able to offer an explanation for cell differentiation 
in vulval development, it cannot capture the full non- 
autonomous nature of the developmental process, since it 
does not reproduce developmental trajectories, chreodes, 
in an accurate and fully time-dependent manner. 

Similarly, most of the few other examples of non- 
autonomous models in the biological literature do not 
explicitly consider the effect of parameter changes on 
transient behaviour (e. g. [91-93]). This simplification may 
be justified in many cases and is necessary for any kind 
of rigorous analytical treatment of a model. In many sit- 
uations, however, it fails to capture essential features of 
the system. For instance, a truly accurate analysis of gap 
gene regulatory dynamics would require the inclusion of 
both non-autonomy from the rapidly changing maternal 
morphogen gradients, and transient dynamics, which are 
known to underlie the temporal shifts in domain posi- 
tion. Before we can undertake such an analysis we must 
first build a conceptual toolkit for phase space analysis of 
transient, non-autonomous dynamics. Due to the limited 
amount of analysis possible in such systems, this toolkit 



will need to be developed in a pragmatic and empirical 
manner, using numerical simulation and exploration of a 
simple conceptual model as its basis. 

In the following sections, we present such a simulation- 
based attempt at developing concepts to classify transient, 
non-autonomous behaviours. For this purpose, we use a 
simple two-component model of a genetic toggle switch, 
whose potential landscape can be explicitly visualised. 
We use time series of graphs and animations of systems 
dynamics on this potential to identify mechanisms leading 
to state transitions, and other forms of pattern forma- 
tion. From this, we are able to identify four basic types of 
dynamical mechanism and behaviour — transitions, pur- 
suits and two types of captures of trajectories — that can 
be used to classify and understand dynamical behaviour in 
more complicated and realistic models, such as a full non- 
autonomous version of the gap gene model. While this 
present paper provides the methodological foundation for 
such an analysis, detailed results for a realistic model of 
the gap genes will be presented elsewhere. 

Model and methods 

To develop our methodology for analysing transient 
behaviour in non-autonomous dynamical systems, we use 
a simple toggle switch model (see [22], and references 
therein) with time-dependent parameters. We consider 
two interacting genes X and Y (Figure IB, panel 1). Con- 
centrations of the corresponding protein products are 
labelled x and y. X and Y mutually repress each other, are 
linearly activated by external signals and can auto-activate 
themselves (Figure IB, panel 1). Protein products decay 
linearly dependent on their concentration. The mathe- 
matical formulation of our toggle switch model is thus 
given by 



dx 
dt 

dy 
dt 



0i X + 



a*+x* 



Otv + 



— 1 



d 4 +x 4 



hvX 



X y y 



(1) 



where parameters a x and a y represent the external acti- 
vation on genes X and Y respectively. Sigmoid functions 
with Hill coefficients of 4 are used to represent auto- 
activation and mutual repression, where parameters a and 
c determine auto-activation thresholds, while b and d 
determine thresholds for mutual repression. Protein decay 
rates are represented by parameters X x and k y . 

The toggle switch model (1) exhibits different dynam- 
ical regimes depending on the values of its parameters 
(Figure 2A-C). Its name derives from the fact that it can 
exhibit bistability over a wide range of parameters. When 
in this bistable region of parameter space, the underly- 
ing phase portrait has two attracting states and one saddle 
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Saddle-node bifurcations 



Pitchfork bifurcation 
(subcritical) 



Figure 2 Dynamical regimes of the toggle switch model. The toggle switch model can exhibit three different dynamical regimes depending on 
parameter values. (A) In the monostable regime, the phase portrait has one attractor point only (represented by the blue dot on the quasi-potential 
landscape). At this attractor, both products of X and Y are present at low concentrations. (B) In the bistable regime, which gives the toggle switch its 
name, there are two attractor points (shown in different shades of blue) and one saddle (red) on a separatrix (black line), which separates the two 
basins of attraction. The attractors correspond to highx, lowy (dark blue), or lowx, high y (light blue). The two factors never coexist when equilibrium 
is reached in this regime. (C) In the tristable regime, both bistable switch attractors and the steady state at low co-existing concentrations are 
present (shown in different shades of blue). In addition, there are two separatrices with associated saddle points (red). These regimes convert into 
each other as follows (double-headed black arrows indicate reversibility of bifurcations): the monostable attractor is converted into two bistable 
attractors and a saddle point through a supercritical pitchfork bifurcation; the saddle in the bistable regime is converted into an attractor and two 
additional saddles in the tristable regime through a subcritical pitchfork bifurcation; the bistable attractors and their saddles collide and annihilate in 
two simultaneous saddle-node (or fold) bifurcations to turn the tristable regime into a monostable one. Graph axes as in Figure 1 B, Panel 4. 



point (Figure 2B). All phase portraits associated with 
parameters in the bistable range are topologically equiva- 
lent to each other, meaning that they can be mapped onto 
each other by a continuous deformation of phase space 
called a homeomorphism [34]. 

The toggle switch model has two other dynamical 
regimes: monostable and tristable. Phase portraits associ- 
ated with parameters in the monostable range have only 
one attractor point (Figure 2A), while those in the tristable 
range have three attractor states and two saddle points 
(Figure 2C). Again, phase portraits within each regime 
are topologically equivalent to each other. While phase 
space can be geometrically deformed within each regime 
(through movements of attractors or separatrices), its 
topology only changes when one regime transitions into 
another through different types of bifurcations [32,34] 
(Figure 2). The transition from monostable to bistable is 
known to be governed by a supracritical pitchfork, the 
transition from bistable to tristable involves a subcritical 



pitchfork bifurcation, and the transition from tristable to 
monostable takes place through two simultaneous saddle- 
node bifurcations involving the two attractors labelled in 
darker blue in Figure 2C. 

Definition of the potential landscape 

Potential landscapes can only be calculated explicitly for 
the class of dynamical systems called gradient systems 
[32] . A two-variable gradient system is a dynamical system 



dx 
dt 



= f(pc,y) 



-dt = g(x ' y) 



(2) 



which satisfies the following relationship between partial 
derivatives 



f y (x,y) =gx(x,y). 



(3) 
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For gradient systems, it is possible to calculate a closed- 
form (explicit) potential function, V(x,y) such that 



Vr 



—dx 
dt 
—dy 



Vy - — 

y dt 



(4) 



The local minima on the two-dimensional potential sur- 
face given by V(x,y) correspond mathematically to the 
steady states of the system in (2) since, if (x*, y*) is such 
that 



(5) 



then 



dx 
dt 
dy 
dt 



(x*,f) = 0 
(**,/) = 0, 



(6) 



and, therefore, (x*, y*) is a steady state of (2). 



Calculating quasi-potential landscapes 

Condition (3) will not always be met. In particular, dynam- 
ical systems representing gene interaction networks are 
not in general gradient systems, and therefore an associ- 
ated potential function and landscape may not exist. In 
such cases, we can still take advantage of the visualisa- 
tion power of potential landscapes by approximating the 
true potential using a numerical method. The numeri- 
cal approximation method we adopt for our study was 
developed by Bhattacharya and colleagues [38] using a 
toggle switch model very similar to the one used here. This 
allows us to calculate a quasi-potential landscape for any 
specific set of fixed parameter values. 

The quasi-potential, which we denote by V qj is defined 
to decrease along all trajectories of the dynamical system 
as they progress on the phase portrait over time 

8V a 8V a —dx —dy , x 

AV q = -^Ax+ ^Ay= — Ax + ~^Ay. (7) 
8x 8y dt dt 

Ax and Ay are defined as small-enough increments 
along the trajectory such that dx/dt and dy/dt can be con- 
sidered constant in the interval + Ax), (y,y+ Ay)]. 
In addition, A* = ^At and Ay = %At, where At is the 
time increment. Substituting into equation 7, we obtain 



AV a = 



—dx / dx ^ 
dt \ dt 



(dx . 
Tt 1 4 



+ 



(—At 
dt \dt 



fdxY 
\di) 



(8) 



At 



AV q has been formulated in such a way that, for pos- 
itive time increments At, AV q is always negative along 
the unfolding trajectory and is, in effect, a Lyapunov func- 
tion of the two-gene dynamical system [32] . This ensures 
that trajectories will always "roll" downhill on the quasi- 
potential surface. Just as in the case of closed-form poten- 
tial, the steady states of the system (x*,y*) correspond 
to the local minima on the quasi-potential surface since 
AV q {x*,y*) = 0. 

We apply the numerical approximation method 
described above to trajectories with various initial points 
on the x-y plane. This yields a sampled collection of tra- 
jectories with quasi-potential values associated to every 
one of their points. Next, we apply the following two 
assumptions, in order to construct a continuous quasi- 
potential surface from this sample of discrete trajec- 
tories [38]: 

1. Two trajectories with different initial conditions that 
converge to the same steady state must also converge 
to the same final quasi-potential level (normalisation 
within basins of attraction). 

2. Two adjacent trajectories that converge to different 
steady states will be taken to start from the same 
initial quasi-potential level (normalisation between 
basins of attraction). 

Finally, interpolation of all the normalised trajecto- 
ries results in a continuous quasi-potential landscape. 
Bhattacharya et al. [38] validated this approach by demon- 
strating that the quasi-potential values of the steady states 
were inversely correlated with their probability of occur- 
rence using a stochastic version of the toggle switch 
dynamical system. 

Approximating non-autonomous trajectories 

As we have argued in the Background Section, we cannot 
generally assume that parameter values remain constant 
over time when modelling biological processes. We take 
a step-wise approximation approach to the change in 
parameter values to address this problem (Figure 3). We 
chose a time increment (step size) as small as possible. 
Parameter values are kept constant for the duration of 
each time step. As a consequence, the associated phase 
portrait will also remain constant during this time inter- 
val, and is visualised for each step by calculating a quasi- 
potential landscape as described in the previous section 
(Figure 3C, top row). 

The smaller the time increments considered, the bet- 
ter we are able to approximate continuous changes in 
parameter values, and the consequent changes to the asso- 
ciated phase portrait and quasi-potential landscape. Such 
accurate approximation allows us to faithfully reproduce 
non-autonomous trajectories produced by models with 
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time 



Figure 3 Numerical approximation of non-autonomous trajectories. (A) Toggle switch network. Red arrows representing auto-activation 
indicate time-dependence of threshold parameters a x and a y (see equation 1 ). (B) Values of auto-activation thresholds a x and a y are altered 
simultaneously and linearly over time. The graph shows the step-wise approximation of a continuous change, in this case, an increase in parameter 
values. Step size is taken as small as computational efficiency allows. (C) During every time step, parameters can be considered constant, and the 
phase portrait and (quasi-)potential landscape are calculated for the current set of parameter values. Trajectories are then integrated over the 
duration of the time step using the previous end point as the current initial condition. The result is mapped onto the potential surface. The four 
panels in (C) show examples of potential landscapes (upper panels) calculated based on sets of parameter values at time points indicated by 
dashed arrows from (B). Important events altering the geometry of the trajectory are indicated. Lower panels show the corresponding 
instantaneous phase portraits with the integrated progression of the trajectory across time steps. See Model and methods for details. 



continuously time- variable parameters. This is done by 
integrating trajectories using constant parameters dur- 
ing each time step, and then using the resulting end 
position in phase space as the initial condition for the 
next time step. The resulting integrated trajectories can 
then be visualised by mapping them from the underlying 
phase plane onto the associated quasi-potential landscape 
as described above. This allows us to track and analyse 
in detail how changes in the phase portrait and quasi- 
potential landscape shape the trajectories as the values of 
the parameters are changing. 



Results and discussion 

Using the toggle switch defined in (1) as a concep- 
tual model, we have set out to identify and catalogue 
mechanisms that affect the shape of trajectories in non- 
autonomous systems. In contrast to previous studies 
(see, for example, [19,20]), we do not assume asymp- 
totic behaviour of the system, but focus on transient 
dynamics — that is, the entire trajectory from initial to 
steady state. 

In this context, we use the term 'mechanism' in an 
unusual, but precisely defined way. We do not mean it to 
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imply any specific biochemical interactions — this is how 
most experimental biologists would use the term. Nor 
does our notion of 'mechanism' depend on any specific 
regulatory structure, such as network motifs (reviewed in 
[94]), or the different stripe-forming mechanisms studied 
in [95]. Instead, we use the term mechanism' in a broader, 
more abstract sense: a mechanism is a causal explana- 
tion of dynamics in terms of the dependence of the flow 
on parameter changes, and how these affect the trajec- 
tories of the system. We will illustrate this very abstract 
definition using a number of concrete examples below. 

Changing parameter values in non-autonomous sys- 
tems can change the underlying phase portrait in two 
main ways: it can alter the geometry or the topology of 
phase space (or both) [32,34]. On one hand, a change 
in geometry means that the phase portrait is still com- 
posed of the same basic elements — the same number, kind 
and relative placement of attractors and their associated 
basins, saddle points, and separatrices — but that the exact 
position, shape, and/or size of these elements has been 
altered. On the other hand, a topological alteration of 
phase space occurs when a change in parameter values 
cause the number or stability of steady states to change. 
For example, small changes in a parameter value can cause 
a pair of attractor and saddle points to appear or disappear 
through a saddle-node or fold bifurcation, or an attractor 
can be turned into a saddle point through a transcriti- 
cal bifurcation. Topological changes to the phase portrait 
are therefore associated with bifurcations as studied in the 
context of autonomous dynamical systems [32-34] . 

In what follows, we examine how both geometrical and 
topological changes in phase space affect transient tra- 
jectories in non-autonomous systems. For this purpose, 
we use a numerical simulation approach with three- 
dimensional visualisation of the quasi-potential landscape 
associated with the phase portrait of our toggle switch 
model (see the Model and methods Section, for details). 
We identify four different basic types of non-autonomous 
mechanisms affecting dynamical behaviour: transitions, 
pursuits and two kinds of captures. 

Transitions 

Out of the four trajectory-shaping mechanisms that we 
have identified in this study, transitions are the most 
familiar. They have been found and described in previ- 
ous approaches to the study of non-autonomous biolog- 
ical systems [78-80,91,92]. Transitions in our context are 
equivalent to Thorns catastrophes' [8], and to what Schef- 
fer has termed critical transitions', defined as a shift of 
the system from one attractor to another when it passes a 
given critical point in parameter space [92]. While Thorn's 
and Scheffer's analyses focus on the resulting steady state 
behaviour of the system, we also consider the transient 
dynamics during a transition. 



An example of a transition in the toggle switch model — 
occurring from the bistable via the tristable to the monos- 
table regime — is shown in Figure 4 (see also Additional 
file 1, Supporting Movie SI). This transition is driven by 
an increase in the value of the auto-activation thresholds 
(a and c in equation 1). Changes in auto-activation have 
been previously used to explain the dynamics of stem 
cell differentiation [20,21,47,48], and patterning by lateral 
inhibition [25]. 

In this case, we consider that the system is already at 
steady state at the outset, in its bistable regime with high x 
and low y (Figure 4A). In other words, we assume that the 
system has had time to converge, or that the initial con- 
dition coincides with the neighbourhood of this attractor. 
During the transition, two bifurcation events take place— 
a subcritical pitchfork bifurcation (Figure 4A to B), and 
two simultaneous saddle-node bifurcations (Figure 4D 
to E) — which make the initial steady state and its asso- 
ciated basin of attraction disappear while creating a new 
attractor at which both factors X and Y are present at low 
concentrations (Figure 4B to E). 

Examining transient dynamics during the transition 
reveals important details that shape the trajectory and 
hence the behaviour of the system. First, we note that 
the trajectory initially at steady state does not necessar- 
ily remain anchored to its attractor (Figure 4C and D). 
As the change in parameter values causes the attractor 
state to move, the trajectory's current state falls behind 
and reacts by travelling towards the moving attractor (see 
Pursuit below). Since the flow rate along the trajectory is 
smaller than the velocity of attractor movement, the sys- 
tem is not able to catch up with the moving steady state. 
Hence, it temporally reverts from asymptotic to tran- 
sient behaviour. Obviously, such a reversal can never be 
observed if we only focus on steady state behaviour. 

Second, we observe a delay between the subcritical 
pitchfork bifurcation creating the new attractor state 
(shown in light blue in Figure 4B) and the system switch- 
ing into that basin of attraction (Figure 4E). This delay 
effect was already noted and described by Thorn in his 
analysis [8]. The change from one basin to another (indi- 
cated by a change from dark to light blue on the trajec- 
tory shown in Figure 4E) only occurs once the tristable 
system undergoes two simultaneous saddle-node bifurca- 
tions, where the two saddles on the separatrices collide 
against the two outer attractors and annihilate each other 
(darker blue attractors in Figure 4D and E). A monostable 
system results from these bifurcations. The trajectory sud- 
denly finds itself in a different basin, and eventually it will 
converge to this new attractor state (Figure 4E). 

In this particular example, the shape of the observed 
transient trajectory does not differ much from its equiva- 
lent in an autonomous monostable system. This is because 
for most of the time the trajectory is in the basin of 
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Figure 4 Transition. A transition indicates the switch of a non-autonomous system from one attractor state to another. Upper panels show 
(quasi-)potential surfaces, lower panels phase portraits as in Figure 3C The progress of time is shown through increasingly dark shading, and by the 
arrow at the bottom of the figure. (A) The system starts off in the bistable regime. The trajectory's initial conditions coincide with the attractor at 
high x, lowy (dark blue). The trajectory is therefore at steady state at the outset. (B) Changes in auto-activation thresholds a and c (equation 1 ) over 
time cause the system to undergo a subcritical pitchfork bifurcation and enter the tristable regime (see also Figure 2). (C, D) The trajectory does not 
switch attractors immediately after the bifurcation occurs. However, it does not remain anchored to its current attractor either. Instead, it is left 
behind by the moving attractor, which it starts to pursue (see also Figures 5 and 6). (E) The system enters the monostable regime as the two 
bistable attractors disappear via two simultaneous saddle-node (or fold) bifurcations. The trajectory suddenly finds itself in an alternative basin of 
attraction, and eventually converges to the new, monostable attractor with low x and y . This change in basins of attraction is represented by a 
change in colour of the trajectory shown in (E). See also Additional file 1, Supporting Movie SI . 



attraction of the attractor at low X and Y concentrations 
(light blue part of the trajectory in Figure 4E). Nev- 
ertheless, most developmental processes require tightly 
controlled timing and therefore the observed temporary 
reversal to transient behaviour and consequent delay in 
attractor switch can be significant. These aspects are 
therefore not negligible if we want to achieve a full under- 
standing of the dynamics of the system. Steady state 
analyses, where more or less instant convergence to the 
new attractor state is assumed, will be limited in this 
regard since they are not able to address the question of 
developmental timing. 

In addition to its effect on timing, the delay in switch- 
ing basins of attraction also provides an explanation for 
the irreversibility of developmental pathways. This phe- 
nomenon has been studied in detail in the context of 
stem cell differentiation. Wang and colleagues [19,20] use 
a non-autonomous toggle switch model similar to ours 
to study the transition from monostable (stem cell) to 
bistable behaviour (differentiated cells). The authors sim- 
ulate this transition in both directions (from stem cell to 
differentiated state and back), and show that two different 
delay effects occur depending on the direction of the pro- 
cess. In other words, forward and backward pathways are 
very different, thus explaining why the system is unable 



to retrace its original differentiation pathway when the 
process is reversed. 

Pursuit 

Transient phenomena need to be taken into account to 
be able to fully understand the dynamical repertoire of a 
system. This becomes especially relevant when consider- 
ing what we call pursuit mechanisms. We now consider 
what happens if the external activation on genes X and 
Y in our toggle switch model (a x and a y in 1) are altered 
over time. In biological terms, this could represent the 
systems response to a changing external signal. This par- 
ticular parameter change affects the geometry but not 
the topology of the phase portrait. The system remains 
in the bistable dynamical regime throughout the whole 
parameter range that we explored. 

Let us first consider an increase in activation strength. 
As external production rates increase over time, the 
attractors move away from the origin along the direc- 
tion of the x and y axes respectively (Figure 5, see also 
Additional file 2, Supporting Movie S2). This outward 
movement of the attractors leaves the position of the 
separatrix unchanged throughout the simulation, keeping 
the location and area of the basins of attraction constant 
over time (Figure 5). Since the separatrix does not move, 
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Figure 5 Pursuit stabilising the direction of a trajectory. Pursuit behaviour results from the movement of an attractor. In this case, both attractor 
and trajectory progress in the same direction, which is therefore stabilised even before the system enters any asymptotic regime. Upper panels 
show (quasi-)potential surfaces, lower panels phase portraits as in Figure 3C The progress of time is shown through increasingly dark shading, and 
by the arrow at the bottom of the figure. (A) Pursuit behaviour shows similarities to autonomous dynamics, since the location of the initial condition 
determines the attractor towards which the system will converge. In this example, parameter changes do not alter the position of the separatrix. 
(B-D) However, the approach of the system towards the attractor is very different than in the autonomous case, since both attractors move away 
from the origin as activation strength (represented by a x and a y in equation 1 ) is increased over time. This leads to an enduring pursuit by the 
trajectory of its moving attractor target. See also Additional file 2, Supporting Movie S2. 



trajectories are not able to change basins. Therefore, the 
initial conditions determine which basin of attraction the 
system will remain in, as they do in autonomous systems. 
However, the approach to the attractor is very different 
in this non-autonomous case, since trajectories are drawn 
towards a moving target. It is this continued pursuit that 
plays the main role in shaping them (Figure 5). Even 
though its associated attractor keeps moving over time, 
this sort of pursuit shows how the general direction of a 
trajectory can be stably maintained in a non-autonomous 
system. 

Two specific pursuit scenarios can be distinguished. In 
the first, the trajectory eventually reaches the moving 
attractor. This occurs if the change in parameter values 
increases the flux around the attractor faster than it moves 
its position, or if further change in parameter values no 
longer alters the position of an attractor— in other words, 
the position of an attractor itself converges to a given 
location over time. Under these conditions, the trajecto- 
ries of a non-autonomous system can show asymptotic 
behaviour and come to rest at steady state. However, it 
is not guaranteed that this should take place. There are 
many imaginable scenarios in which attractors will keep 
moving as parameters change and trajectories governed 



by pursuit mechanisms may never come to rest. It is prob- 
able that many developmental processes— at least to a 
certain degree— work in this regime, as external condi- 
tions and the cellular environment (signalling inputs or 
tissue context) constantly keep changing over time. 

Just as in the case of transitions, it is obvious that pur- 
suit mechanisms have a great impact on the timing of 
developmental dynamics. They can delay or even pre- 
vent the system from reaching steady state. In addition, 
attractor movements can also drastically influence and 
alter the shape of a trajectory. This becomes evident if 
we consider an alternative scenario, in which activating 
inputs decrease over time (Figure 6, see also Additional 
file 3, Supporting Movie S3). In this case, the attractors 
move towards the origin over time. This causes the attrac- 
tor to overtake' the trajectory at a given point in time 
(Figure 6B), which induces a rather drastic change in the 
trajectory's direction: initially it is travelling toward high x 
(Figure 6A), but later comes to approach the origin (low#) 
while pursuing the moving attractor (Figure 6B-D). 

In summary, we have illustrated how pursuit mech- 
anisms in non-autonomous systems can either lead to 
the stable maintenance of a transient trajectory (as in 
Figure 5), or to a drastic change in its direction (as in 
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Figure 6 Pursuit altering the direction of a trajectory. Pursuit behaviour results from the movement of an attractor. In this case, the direction of 
the trajectory is altered since the attractor moves against the flow. Upper panels show (quasi-)potential surfaces, lower panels phase portraits as in 
Figure 3C The progress of time is shown through increasingly dark shading, and by the arrow at the bottom of the figure. (A-D) Panels show the 
movement of the attractors towards the origin as activation strength (represented by a x and a y in equation 1) is decreased overtime. As the 
attractor 'overtakes' the trajectory (between panels B and C), the angle of the flow changes drastically, leading to a reversal in the direction of the 
trajectory (clearly visible in D). See also Additional file 3, Supporting Movie S3. 



Figure 6) depending on the geometric arrangement of 
the attractor with regard to the converging trajectory. 
We have argued that pursuit mechanisms are likely to be 
almost ubiquitous in systems that are exposed to a chang- 
ing environment or variable tissue context (e. g. signalling 
or external input by morphogen gradients). One specific 
example for pattern formation by pursuit are the tran- 
sient gap domain shifts examined by Manu and colleagues 
[35,89]. In this particular case, spatial shifts in domain 
position are caused by the timed up- and later down- 
regulation of a specific sequence of gap genes in nuclei 
located within the posterior region of the Drosophila 
blastoderm embryo (see [35], for details). Basically, each 
nucleus goes through a stereotypical series of expression 
pulses of different gap genes. Such pulses can easily be 
explained by the sort of trajectory observed in Figure 6, 
where an initial expression trend (upregulation of X) is 
later smoothly reversed (towards downregulation of the 
same gene). 

Capture 

In autonomous dynamical systems, trajectories never 
cross a separatrix [32]. Depending on the position of its 
initial conditions, a trajectory will find itself in a particu- 
lar basin of attraction, that is, in the dynamical regime of 



a particular attractor. The trajectory will always remain in 
this basin and, given enough time, will eventually converge 
to the attractor. 

The same is not true for non-autonomous systems. We 
know that the phase portrait is changing as the parame- 
ters change. This means that the position and number of 
steady states need not remain constant over time. We have 
seen how saddles and attractors can be created and anni- 
hilated by bifurcations, and how they can move through 
phase space in such systems. In particular, when saddle 
points change their position, they cause their associated 
separatrices to shift with them. Separatrices mark the 
boundaries between basins of attraction, and their move- 
ment will cause the size, shape and/or relative placement 
of these basins with regard to each other to change. 

As trajectories progress on the changing phase portrait, 
it is possible for a moving separatrix to overtake them. We 
will call this event a capture! When a capture occurs, the 
affected trajectory is deviated, and often exhibits a sud- 
den change in direction, since it is now attracted towards 
a different attractor at a different position of the phase 
plane than it was before. Captures will happen when the 
flow along the trajectory is smaller than the rate of change 
in separatrix position. The separatrix is recruiting points 
from one basin of attraction into another at a faster rate 
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than the trajectory is travelling away from it. Therefore, 
the velocities of the trajectory and the separatrix relative 
to each other are determining whether a capture event is 
going to take place or not. 

We identified two main ways by which capture events 
are brought about in the non-autonomous toggle switch 
model. In both cases, the mechanism shaping the trajec- 
tory is essentially the same: the trajectory suddenly finds 
itself converging towards a different attractor as a result 
of an abrupt change to the flow at its current position. 
The main difference between the two capture mechanisms 
resides in how the movement of the separatrix is caused. 

In the first situation, the trajectory gets captured after 
a bifurcation event has lead to the creation of a new 
attractor state. This increases the number of attractor 
basins and, in this way, introduces new separatrices into 
the phase portrait. We simulate this situation using the 
bistable-to-tristable transition caused by an increase in 
the auto-activation threshold as described above (Figure 7, 



see also Additional file 4, Supporting Movie S4). In this 
example, a subcritical pitchfork bifurcation creates a new 
attractor and two associated saddles from a pre-existing 
saddle point (Figure 7A,B). This results in a change in 
phase space topology. What used to be a single separatrix 
now opens up', giving rise to two different forked sep- 
aratrices (Figure 7B,C). Further parameter changes then 
cause the new separatrices to move outward through 
phase space, catching up with, and overtaking, trajecto- 
ries as they recruit points into the newly created and 
expanding basin of attraction (Figure 7B-D). 

As with the transition mechanism, we notice a some- 
times considerable time delay between the bifurcation and 
the capture event. The extent of this delay depends on 
how close a trajectory is to the new attractor at the time of 
bifurcation, and on the rate of change in separatrix posi- 
tion due to further parameter change. In our example, 
the delay effect can be clearly seen, as the trajectory first 
veers to the left, towards its original attractor, before being 
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Figure 7 Capture due to a change in the topology of the phase portrait. A capture results from a trajectory being recruited into a new basin of 
attraction due to the movement of a separatrix. In this example, the relevant separatrix is created and caused to move by a preceding bifurcation 
event, which leads to the appearance of a new attractor state, resulting in a change of phase space topology. Upper panels show (quasi-)potential 
surfaces, lower panels phase portraits as in Figure 3C. The progress of time is shown through increasingly dark shading, and by the arrow at the 
bottom of the figure. (A) The system starts off in the bistable regime and the initial conditions place the trajectory in the basin of the high x, lowy 
attractor (dark blue). (B) Changes in the values of the auto-activation thresholds {a and c, see equation 1) cause the system to undergo a subcritical 
pitchfork bifurcation and enter the tristable regime (see also Figure 2). At the time of the bifurcation, the trajectory is still attracted towards the dark 
blue attractor. (C) As auto-activation thresholds are further increased, the two separatrices surrounding the new attractor state (shown in light blue) 
separate from each other, enlarging the corresponding basin of attraction. A capture event takes place as the separatrix 'overtakes' the trajectory, 
recruiting it into the new basin of attraction. The system will now converge towards the light blue attractor. This change in basins of attraction is 
represented by the colour coding of the trajectory on the phase portrait. (D) As auto-activation thresholds are further increased, the system will 
transition from the tristable into the monostable regime (see also Figures 2 and 4). This causes the dark blue attractors and their basins to disappear 
altogether, but does not influence the direction of the trajectory anymore, which will eventually converge to the light blue attractor at low x and y. 
See also Additional file 4, Supporting Movie S4. 



Verd et al. BMC Systems Biology 201 4, 8:43 
http://www.biomedcentral.eom/1752-0509/8/43 



Page 15 of 19 



steered back towards the center of the phase portrait, 
where the new attractor lies. This introduces a marked 
bend into the trajectory (Figure 7D), which affects both 
time to convergence and the identity of transient states. 

In the second scenario, we consider a situation where 
the capture event is not preceded by a bifurcation. In this 
case, a pre-existing separatrix is moving through phase 
space due to parameter changes. This movement of the 
separatrix reconfigures the geometry of the basins of 
attraction without changing the topology of the phase por- 
trait. We simulate this event in our toggle switch model 
by introducing asymmetric changes in the values of the 
thresholds determining the two mutually repressing reg- 
ulatory interactions (i.e. b ^ d in 1) (Figure 8, see also 
Additional file 5, Supporting Movie S5). This shifts the 
position of the separatrix in the bistable regime towards 
one of the two attractors without creating or annihilating 
any steady states. 



As in the case of capture with preceding bifurcation, the 
separatrix moves and recruits points from one basin into 
the other, thereby shrinking one basin of attraction while 
enlarging the other (Figure 8A-D). If the moving sepa- 
ratrix encounters a trajectory at one of the points that 
are being recruited into a different basin of attraction, a 
capture event will take place (Figure 8B). The affected tra- 
jectory will change direction as it now converges towards 
a different attractor than before the capture (Figure 8C,D). 

Capture mechanisms illustrate how dramatic changes in 
direction and velocity of trajectories do not always require 
changes to the topology of the underlying phase portraits. 
In other words, they do not necessarily have to be associ- 
ated with bifurcations, and even if they are, the observed 
sudden change in direction does not usually coincide with 
the bifurcation. Instead, it can be significantly delayed. 
This kind of phenomenon cannot be observed if only 
asymptotic behaviours are considered. In that case, only 



A Potential Landscape 




time 



Figure 8 Capture due to a change in the geometry of the phase portrait. A capture results from a trajectory being recruited into a new basin of 
attraction due to the movement of a separatrix. In this example, the relevant separatrix is caused to move by a break in the symmetry of the 
repressive interactions between X and Y. No new attractor states or separatrices are created. Therefore, the change in landscape topography is 
purely geometrical and does not affect the topology of phase space. Upper panels show (quasi-)potential surfaces, lower panels phase portraits as 
in Figure 3C. The progress of time is shown through increasingly dark shading, and by the arrow at the bottom of the figure. (A) The system starts 
off in the bistable regime and the initial conditions place the trajectory in the basin of the low x, high y attractor (light blue). (B) Changing the 
threshold of one of the repressive interactions only {b in equation 1) shrinks the basin of attraction of the lowx, high y attractor (light blue) causing 
the separatrix to move towards the upper left corner of the phase portrait. (C) The shifting separatrix catches up with the converging trajectory, 
recruiting it into the basin of the high x, lowy attractor (dark blue). A capture event has taken place without any preceding bifurcation. This change 
in basins of attraction is represented by the colour coding of the trajectory on the phase portrait. (D) After the capture, the system converges to its 
new attractor (dark blue). Interestingly, the geometry of this capture has made the trajectory loop over itself. Such self-crossing trajectories are never 
observed in autonomous dynamical systems. The subsequent saddle-node bifurcation by which the saddle (red) and the light blue attractor 
annihilate each other does not affect the trajectory any further. See also Additional file 5, Supporting Movie S5. 
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transition-like events can occur since captures and pur- 
suits crucially depend on transient dynamics. 

There are several examples which illustrate how this 
might be important in practice. Corson and Siggia [90] 
considered an explicitly non-autonomous mechanism for 
determining vulval cell fates in their model. In their study, 
shifts in separatrices and captures happen instantaneously 
at the outset of the simulation. In other words, the initial 
condition of cells that are captured completely determines 
which basin of attraction they will find themselves in for 
the rest of the simulation. While this does not affect the 
analysis of the resulting differentiated (steady) states — 
as evidenced by the tight match between the model and 
the measured proportions of differentiated cells [90] — it 
could limit the models ability to reproduce the transient 
dynamics of cellular differentiation. 

This becomes crucial when analysing more dynamic 
examples of gene expression, such as pattern formation 
driven by the gap gene system. In this example, it will be 
essential to study the interplay between changing mater- 
nal gradients and moving gap target domains. This was 
not really possible in the study of Manu and colleagues 
[35] since their model assumed the maternal gradients to 
remain constant. Our classification scheme of transient 
behaviour in non-autonomous systems should be useful to 
analyse such questions in more detail. For instance, cap- 
ture events occurring at different times in different nuclei 
could explain how stable expression boundaries can be 
maintained in the presence of changing regulatory inputs. 
In contrast, moving domain boundaries in the posterior of 
the embryo are more likely to involve some type of pur- 
suit mechanism (see above). Without doubt, the concepts 
developed in this paper will be useful to classify and char- 
acterise the boundary-forming and -shifting mechanisms 
in a fully non-autonomous model of this system. 

Conclusions 

In this paper, we have argued that considering both 
transient dynamics and explicit time-dependence (non- 
autonomy) of dynamical systems is essential for an accu- 
rate and complete quantitative understanding of many 
biological processes. In particular, we have argued that 
these two features are important if we are to fully cap- 
ture the richness and explanatory power of Waddingtons 
metaphor of the epigenetic landscape [3]. Waddingtons 
concepts of developmental chreodes, and their canalised 
nature due to the regulatory process of homeorhesis, 
apply to transient or non-autonomous behaviour, not sta- 
ble steady states of a system. Moreover, Waddingtons 
landscape changes its shape due to genetic or environ- 
mental signals as represented by the pegs and ropes in 
Figure 1A. 

However, it remains difficult to analyse transient 
dynamics in non-autonomous systems. Where rigorous 



analysis is possible, it is mostly limited to special cases 
(e.g., systems with regular external forcing, [82]). No gen- 
eral theory or toolkit exists to deal with more complex 
non-autonomous systems, or the transient phenomena 
they exhibit. 

In this paper, we have taken a first step towards 
overcoming this limitation. We have used a pragmatic 
simulation-based approach to classify different mecha- 
nisms that affect transient trajectories in non- 
autonomous systems. Our study is based on a simple 
conceptual model— a toggle switch with time-dependent 
parameter values — to visualise the changing potential 
landscape which determines the dynamical repertoire 
of a regulatory circuit. We found four distinct mecha- 
nisms that affect trajectories on time- variable potential 
surfaces: transitions, pursuits and two kinds of captures. 
At this point it is impossible to say whether this list is 
complete. Other types of behaviour may be possible in 
higher-dimensional phase spaces. 

Despite its potential incompleteness, we believe that 
our preliminary classification scheme is useful, since it 
makes transient dynamics in non-autonomous systems 
more tractable in practice. The advantages of our frame- 
work are twofold: first, it allows us to describe and 
categorise patterning mechanisms. Each one of the mech- 
anisms described above has its own characteristics, both 
in terms of their effect on the short-term dynamics, as well 
as the longer-term evolutionary potential of the system. 
We have described specific examples of such character- 
istics in the Results Section. Second, it can be used as 
a diagnostic framework for the analysis of phase space, 
even in the case of higher-dimensional systems where the 
potential surface can no longer be visualised explicitly. 
It shows us, for example, that we need not necessarily 
find a bifurcation associated with every observed critical 
transition of a system. Instead, drastic changes in systems 
behaviour can also be caused by pursuit or geometrical 
capture events as described in the Results Section. 

Our framework is widely applicable to the analysis 
of real- wo rid regulatory networks. In fact, it has been 
developed in the context of our analysis of a fully non- 
autonomous version of the gap gene model by Manu et al 
[35]. Our classification scheme has allowed us to cate- 
gorise different mechanisms for spatial pattern formation 
in this system, and to identify those features of phase space 
responsible for the spatio-temporal dynamics of gene 
expression. A full analysis of this system will be published 
elsewhere. We hope that our work will inspire further such 
analyses in systems where our framework could be useful, 
and where the required detailed and data-driven models 
are currently available. Examples such as vulval determi- 
nation in C. elegans or neural patterning in vertebrates 
come to mind. We have no doubt that— considering the 
current rate of progress in systems biology— this list will 
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be extended considerably in the very near future, leading 
to a more general investigation of transient dynamics in 
non-autonomous biological systems. 

Additional files 



Additional file 1 : Supporting Movie SI -Transition. A transition 
indicates the switch of a non-autonomous system from one attractor state 
to another. This movie shows the simulation that Figure 4 is based on. It 
displays the changing potential landscape as the system transitions from 
the bistable to the tristable and then the monostable regime (see Figure 2) 
as auto-activation thresholds are increased. The corresponding phase 
portrait is shown in the inset on the upper right. Parameter values are 
displayed in the panel on the lower right. The current state of the system is 
shown as a black ball. Attractors are marked by blue dots, saddle nodes by 
red dots with associated separatrices shown as black lines. The history of 
the trajectory is shown on the phase portrait, colour coded according to 
which attractor the trajectory is converging towards at every individual 
step of the simulation. See main text and Figure 4 for details. 

Additional file 2: Supporting Movie S2-Pursuit stabilising the 
direction of a trajectory. Pursuit behaviour results from the movement of 
an attractor. In this case, both attractor and trajectory progress in the same 
direction, which is therefore stabilised even before the system enters any 
asymptotic regime. This movie shows the simulation that Figure 5 is based 
on. It displays the changing potential landscape as activation strength is 
increased. The corresponding phase portrait is shown in the inset on the 
upper right. Parameter values are displayed in the panel on the lower right. 
The current state of the system is shown as a black ball. Attractors are 
marked by blue dots, the saddle node by a red dot with its associated 
separatrix shown as a black line. The history of the trajectory is shown on 
the phase portrait. See main text and Figure 5 for details. 

Additional file 3: Supporting Movie 3-Pursuit altering the direction 
of a trajectory. Pursuit behaviour results from the movement of an 
attractor. In this case, the direction of the trajectory is altered since the 
attractor moves against the flow. This movie shows the simulation that 
Figure 6 is based on. It displays the changing potential landscape as 
activation strength is decreased. The corresponding phase portrait is 
shown in the inset on the upper right. Parameter values are displayed in the 
panel on the lower right. The current state of the system is shown as a black 
ball. Attractors are marked by blue dots, the saddle node by a red dot with 
its associated separatrix shown as a black line. The history of the trajectory 
is shown on the phase portrait. See main text and Figure 6 for details. 

Additional file 4: Supporting Movie 4-Capture due to a change in the 
topology of the phase portrait. A capture results from a trajectory being 
recruited into a new basin of attraction due to the movement of a 
separatrix. In this example, the relevant separatrix is created and caused to 
move by a preceding bifurcation event, which leads to the appearance of a 
new attractor state, resulting in a change of phase space topology. This 
movie shows the simulation that Figure 7 is based on. It displays the 
changing potential landscape as the system transitions from the bistable 
to the tristable and then the monostable regime (see Figure 2) as 
auto-activation thresholds are increased. The corresponding phase portrait 
is shown in the inset on the upper right. Parameter values are displayed in 
the panel on the lower right. The current state of the system is shown as a 
black ball. Attractors are marked by blue dots, saddle nodes by red dots 
with associated separatrices shown as black lines. The history of the 
trajectory is shown on the phase portrait, colour coded according to which 
attractor the trajectory is converging towards at every individual step of 
the simulation. See main text and Figure 7 for details. 

Additional file 5: Supporting Movie 5 — Capture due to a change in 
the geometry of the phase portrait. A capture results from a trajectory 
being recruited into a new basin of attraction due to the movement of a 
separatrix. In this example, the relevant separatrix is caused to move by a 
break in the symmetry of the repressive interactions between X and Y. No 
new attractor states or separatrices are created. Therefore, the change in 
landscape topography is purely geometrical and does not affect the 



topology of phase space. This movie shows the simulation that Figure 8 is 
based on. It displays the changing potential landscape as the threshold for 
repression of X by Y is raised, causing the system to be increasingly 
asymmetrical. The corresponding phase portrait is shown in the inset on 
the upper right. Parameter values are displayed in the panel on the lower 
right. The current state of the system is shown as a black ball. Attractors are 
marked by blue dots, the saddle node by a red dot with its associated 
separatrix shown as a black line. The history of the trajectory is shown on 
the phase portrait, colour coded according to which attractor the 
trajectory is converging towards at every individual step of the simulation. 
See main text and Figure 8 for details. 
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